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1.  Introduction 

Surface  wave  magnitudes  are  an  indispensable  tool  for  discriminating  between  shallow 
earthquakes  and  explosions,  at  least  down  to  body  wave  magnitude  mb  -  4.0.  Seismologists  have 
exhaustively  studied  surface  wave  magnitudes  for  years  (e.g.,  Gutenburg,  1945;  Vanek  et  al., 
1962;  von  Seggem,  1977;  Okal,  1989;  Rezapour  and  Pearce,  1998;  and  many  others).  However, 
operational  methodologies  for  measuring  surface  waves  still  rely  primarily  on  measuring 
unfiltered  dispersed  surface  waves  in  the  time  domain,  primarily  in  the  vicinity  of  20-second 
periods.  A  major  problem  with  this  approach  is  the  significant  effect  of  non-dispersed  Rayleigh 
waves  (Airy  phases)  which  can  occur  at  both  regional  and  teleseismic  distances,  and  can  occur 
with  dominant  periods  much  less  than  20  seconds.  Errors  introduced  by  measuring  Airy  phase 
signals  can  be  greater  than  0.5  magnitude  units,  or  greater  than  a  factor  of  three,  which  is 
unacceptable  for  reliable  network  averaging.  Methods  have  been  developed  (Marshall  and 
Basham,  1972)  which  apply  empirical  corrections  based  on  geological  regions  for  signals  less 
than  20  seconds,  but  this  only  partially  addresses  the  Airy  phase  problem. 

With  digital  processing  now  widely  available,  we  can  optimize  time  domain  processing  and 
make  possible  automated  routine  measurements  at  variable  periods.  This  paper  investigates  the 
use  of  narrow-band  Butterworth  filters  for  measuring  surface  waves,  implemented  as  simple  time 
domain  digital  filters,  and  pays  attention  to  the  effect  of  filtering  on  amplitudes  and  dispersion. 
The  starting  point  is  Herrmann  (1973),  who  examined  in  detail  theoretical  surface  wave  envelope 
functions  derived  from  rectangular  and  Gaussian  frequency  domain  filters,  and  Yacoub  (1983), 
who  applied  this  methodology  to  automating  surface  wave  measurements  between  17-23 
seconds,  using  narrow-band  Gaussian  filters.  I  also  demonstrate  how  to  modify  currently 
accepted  surface  wave  magnitude  formulae  to  be  unbiased  with  respect  to  narrow-band  filtering 
at  variable  periods. 

The  first  portion  of  the  paper  focuses  on  transformations  of  narrow-band  Butterworth  filters 
operating  on  dispersed  waveforms  from  the  frequency  domain  to  the  time  domain,  without 
incorporating  geometric  spreading  and  attenuation  effects.  Asymptotic  formulae  are  evaluated 
for  simple  frequency/time  transformations,  with  corresponding  error  analysis. 

The  second  part  of  the  paper  reviews  surface  wave  magnitudes  from  a  theoretical  point  of  view, 
by  correcting  frequency  domain  amplitudes  for  dispersion,  geometric  spreading,  and  attenuation, 
and  showing  how  this  results  in  currently  accepted  magnitude  corrections.  Using  the  results  of 
the  first  part,  narrow-band  Butterworth  measurements  are  incorporated  into  accepted  magnitude 
corrections,  resulting  in  a  new  formula  which  is  unbiased  with  respect  to  20-second 
measurements,  at  variable  periods  between  5-25  seconds. 
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2.  Surface  Wave  Amplitudes 

A  filtered,  dispersed  propagating  surface  wave  normal  mode  can  be  expressed  as  (Herrmann, 
1973): 


a(t,x)  =  —  J  H(co)A(co)eiico,~ki<0)x)dco 
2iz 


a) 


where  tois  the  angular  frequency,  H(a>)  is  a  band-pass  filter  symmetric  about  center  frequency 
cco,  A(co)  is  the  complex  amplitude  of  the  normal  mode,  and  k(co)  is  the  wavenumber.  A(co)  is 
also  distance  dependent  if  geometric  spreading  and  attenuation  are  considered.  These 
corrections  will  be  ignored  for  now  to  concentrate  on  the  effects  of  dispersion,  and  discussed  in 
detail  in  the  section  on  surface  wave  magnitudes.  To  approximate  the  dispersive  characteristics 
of  the  propagating  wave,  expand  the  phase  as  a  Taylor  series  about  a  center  frequency  OX), 
ignoring  higher  order  terms  beyond  the  quadratic: 

(cot-kx)  =  </>0+fi(m-W0)-a(m-m0)2  (2) 


where  (Herrmann,  1973;  Aid  and  Richards,  1980) 
<f>0=G)0t-k0x,  =  t-  x——\  =  t- 


dm 


Uc 


a- 


x  d2k 


2  dco2 


x  r0  du 


Aitul  dT 


U  is  the  group  velocity  and  T  is  the  period  corresponding  to  the  angular  frequency  ax  Following 
Papoulis  (page  123,  1962)  we  make  the  assumption  that  H(co)  is  a  narrow-band  symmetric  filter 
and  A(co )  is  approximately  constant  across  the  bandwidth  of  H.  Then,  substituting  (2)  into  (1) 
and  making  the  variable  substitution  co+  axo  —>  QX  gives  the  following: 


*((,*)  =  V*" 1  f  (4) 

7t  J 


where  A0  =  A(a>o),  and  Hl(cd)  =  H(co+  m)  is  an  equivalent  low-pass  filter.  The  real  part  of  the 
complex  expression  (4)  is  equivalent  to  (1).  The  envelope  maximum  corresponding  to  the  group 
velocity  p =  0  is  defined  as: 


Mo1 

oo 

f  H L{a>)e~iac°2  dm 

1  n 

J 

—CO 

(5) 


This  expression  shows  that,  for  a  narrow-band  process,  the  maximum  time  domain  amplitude  is 
equivalent  to  the  frequency  domain  amplitude  modulated  by  two  low  pass  filters:  Hi  and 
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e  tao}  .  For  values  of  ©away  from  the  origin,  the  exponential  integrand  will  rapidly  oscillate, 
not  contributing  to  the  integral,  and  thus  act  as  a  low  pass  filter.  As  a  result,  the  time  domain 
amplitude  can  be  controlled  by  either  filter,  depending  on  the  value  of  the  cutoff  frequency  an,  of 
Hi,  or  the  value  of  or  in  the  exponent.  This  can  be  seen  graphically  in  Figure  1 : 


“0.04  -0.02  0  0.02  0.04 


l-0.05  f.  JD.05 

Frequency  (Hz) 

Alpha  Exponent 
'  ‘ '  Low  Pass  Filter  HL 

Figure  1. 

As  the  value  of  a  increases,  the  width  of  the  fundamental  lobe  of  the  exponential  integrand  will 
decrease,  thus  controlling  the  bandwidth  of  the  recorded  amplitude.  If  the  value  of  the 
bandwidth  of  Hi  decreases  below  the  exponential  bandwidth,  it  will  control  the  amplitude  of  the 
recorded  signal.  This  trade-off  can  be  useful  in  modifying  the  dispersion  measured  on  observed 
seismograms,  as  will  be  discussed  below. 

3.  Application  of  Butte rworth  Filters 

As  a  specific  case  of  the  low-pass  filter  Hi,  let  us  use  a  zero  phase  n*  order  Butterworth  band¬ 
pass  filter  of  the  form  (Kanasewich,  1975) 

"<■»>—;  '  s2„  («) 

1+  0 
coc 

\  c  J 

where  2©  is  the  bandwidth  about  COq  ,  defined  at  one-half  the  amplitude  of  H(o).  Then, 
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Hl(co)  = 


1+  -2- 


,  ,,2« 

CO  +  C0c 


Substituting  (7)  into  (4)  gives 


a(t,x)  =  A0e  -  J  — - —  e 

x  J  CO  +  cot 


UPm-ouo1)  d0) 


Expanding  Hl  as  partial  fractions  and  using  relationship  7.4.2  from  Abromowitz  and  Stegun 
(1964),  a(t,x)  can  be  shown  to  equal  (see  Appendix  A): 


«((,*)  =  V**0  (9) 

In  4—  L  J 


where  erfc  is  defined  as  the  complementary  error  function  (l-erf), 


%(±fi)  =  e^ek±i^= 


£k=<»ce 


jt{2k-n-\) 


and  recalling  from  equation  (3)  that 


n  A 

B-t - , 

U0 


±_T^_dU_ 

4 n  ul  dT  o 


As  an  example,  set  (typical  for  continental  crusts) 


X=  6000  km,  Tq=  20  sec,  UQ=  2.9  km/sec,  —  =.02  km/sec 

dT  n 
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Let  the  order  of  the  Butterworth  filter  be  n  =  3,  and  the  comer  of  the  filter  fc  =  .005  Hertz. 
Substituting  these  values  into  equation  (9)  and  plotting  the  real  value  as  a  function  of  /?  gives  the 
seismogram  illustrated  in  Figure  2. 


Filtered  Surface  Wave 


-400  -300  -200  -100  0  100  200  300  400 

t-x/U 
Figure  2. 


Notice  that  the  plot  has  a  maximum  at  the  group  velocity  Uo  =  x/t. 

4.  Asymptotic  Values  of  Filtered  Amplitudes 

At  the  amplitude  maximum  corresponding  to  the  group  velocity  (fi  =  0),  equation  (9)  reduces  to: 

a0  =  A0el</,°  -\£kel0l£k  erfc(Jiaek)  (15 

n 


where 


Sk=CQce  Uk 


7t(2k  -n- 1) 


For  large  (ti  or  a,  the  following  asymptotic  form  for  erfc  can  be  used  (Abromowitz  and  Stegun, 
relation  7.1.23,  1964;  Mathews  and  Walker,  page  80,  1970): 


erfc(u) 
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Letting  u  =  £k ,  and  substituting  (13)  and  (18)  into  (15)  gives,  after  some  algebraic 
manipulation, 


Ao 

jaicc 

71, 1 

\dU 

X 

1 

dT 

(19) 


This  is  precisely  the  value  arrived  at  by  Okal  (1989)  assuming  strongly  dispersed  surface  waves 
(a  large).  Notice  that  equation  (19)  is  independent  of  cq,  thus  independent  of  Hi. 

For  small  values  of  <2  in  the  exponential  integrand,  it  is  expected  that  the  main  lobe  will  be  wide 
(Figure  1),  so  the  amplitude  will  be  controlled  by  the  low-pass  filter  Hi.  Letting  a  0  in 
equation  (15)  results  in 

«0=VA— e°> 

n  *=i 


which  is  controlled  only  by  the  order  of  the  Butterworth  filter  Hl  and  its  comer  frequency  cq. 

For  a  3rd  order  filter  (n  =  3),  and  recalling  equation  (17),  equation  (20)  immediately  reduces  to 

=  <2!> 


Figure  3  shows  the  envelope  of  equation  (15)  as  a  function  of  ct,  along  with  the  asymptotic 
values  (19)  and  (21),  assuming  a  3rd  order  Butterworth  filter  with  a  filter  comer  of 
fc  =  .005  Hertz: 


Alpha 


Butterworth  Filter 
‘ '  Small  Alpha  Asymptote 
— '  Large  Alpha  Asymptote 

Figure  3. 
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Notice  that  the  Butterworth  filtered  surface  wave  shows  almost  constant  amplitude  (minimal 
dispersion)  out  to  a  value  of  a=  800,  where  the  asymptotes  intersect. 

For  any  value  of  a,  cq,  and  n,  this  can  be  generalized  as  follows:  looking  only  at  absolute  values, 
rewrite  (19)  and  (21)  as 


(22) 


A0O)crn 


=  -£cos(0j 

n  *=l 


(23) 


Notice  that  the  exponential  term  in  r„  reduces  to  the  cosine  term  due  to  the  symmetric 
distribution  of  0*  (see  equation  17).  Define  the  asymptotic  intercept  as  the  point  where  the 
equations  are  equal: 


cor  = 


■  1 


not 


(24) 


Assuming  minimal  dispersion  for  values  of  a  less  than  the  intercept  results  in 

or  equivalently, 

co  < — \= 
r^iza 


(25) 


Equation  (25)  determines  the  range  of  cutoff frequencies  for  which  a  Butterworth  filter  can  be 
constructed  to  reduce  the  dispersion  error  to  a  minimum. 

To  justify  the  assumption  of  minimal  dispersion,  the  absolute  error  between  the  Butterworth 
amplitude  and  the  asymptotic  value  at  the  intercept  can  be  determined  as  follows.  Evaluate  and 
rearrange  the  Butterworth  amplitude  (15)  at  the  intercept  point  (24)  as: 

aiNT  ~  ^0  ^/J  (26) 


where 


erfc(Vk), 


It  can  be  seen  by  inspection  that  Z'„  is  only  a  function  of  the  Butterworth  order  n. 


(27) 
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Define  the  logarithmic  error  between  the  Butterworth  amplitude  and  asymptote  as: 

ERR  =  LOG(ASYMPTOTE)  -  LOG(AMPLITUDE) 
Evaluate  at  the  intercept  using  (23)  and  (26)  for: 

ERR  =LOG(A0a>cr„)-LOG{A0a>c  Z„) 


or 


ERR  =  LOG 


'i' 

\Z“  J 


(28) 


Notice  that  the  logarithmic  error  is  only  a  function  of  the  Butterworth  order  n.  Table  1  gives  the 
values  of  ERR  for  the  first  six  Butterworth  orders: 


Table  1. 

Butterworth  Order  vs.  Log  Error 


n 

1 

2 

3 

4 

ERR 

0.191 

0.0672 

0.0360 

0.0244 

For  Butterworth  filters  of  order  3  (recommended)  or  greater,  the  maximum  magnitude  error  will 
be  less  than  or  equal  to  0.0360  for  values  of  cck,  satisfying  the  inequality  in  (25). 

5.  Surface  Wave  Magnitudes 

Surface  wave  magnitudes  are  generally  expressed  as  distance-corrected  logarithmic  amplitude 
measurements  of  surface  waves  in  the  time  domain,  usually  measured  in  the  vicinity  of 
20-second  vertical  component  Rayleigh  waves.  Typically,  the  amplitude  measurements  are 
corrected  for  instrument  response  to  either  zero-to-peak  or  peak-to-peak  amplitude 
measurements,  usually  in  millimicrons  or  nanometers.  Time  domain  formulae  for  surface  wave 
magnitudes  are  generally  derived  from  empirical  measurements  of  surface  wave  amplitudes 
averaged  across  many  events  and  epicentral  distances  (e.g.,  Vanek  et  ai.,  1962;  von  Seggem, 
1977;  Rezapour  and  Pearce,  1998). 

Surface  wave  magnitudes  can  also  be  derived  theoretically  by  correcting  frequency  domain 
amplitudes  for  geometric  spreading  and  attenuation  (Kanamori  and  Stewart,  1976)  and  then 
transforming  to  the  time  domain  (Okal,  1989).  Let 


TttcS 

Ac  =  A*Jre  sin(A)  eUQT 


(29) 


where 
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Ac  = 

Corrected  frequency  domain  amplitude 

A  = 

Frequency  domain  amplitude 

re  = 

Earth’s  radius 

A  = 

Epicentral  distance  in  degrees 

K 

Degree  to  kilometer  conversion  (111.2  km/deg) 

T  = 

Period  of  interest 

Q  = 

Q  factor  measuring  attenuation  at  period  T 

u  = 

Group  velocity  at  period  T 

To  transform  (29)  into  the  time  domain,  recall  the  strong  dispersion  relation  for  surface  wave 
amplitudes  ( 1 9)  (see  also  Okal,  1989): 


a  = 


2U 


A 


(30) 


Solve  (30)  for  A  and  substitute  into  (29)  to  arrive  at  a  theoretical  expression  for  corrected  time 
domain  amplitudes  at  frequency  T: 


, _  -I 

ac  =  aT  ^Jresin(A)eUQT  — 


dU 

dT 


kA 


2  U 


(31) 


Taking  the  base  10  logarithm  of  (31)  results  in  an  expression  for  surface  wave  magnitudes: 


Ms  =  log (aT) + |log(sin(A)) + log(e)  + |log(A)  +  logj 


(U 


+  C 


(32) 


where  constants  are  now  combined  in  C.  The  correction  terms  represent  distance  adjustments  for 
geometric  spreading,  attenuation,  and  dispersion.  The  last  term  is  a  period-dependent  dispersion 
correction.  The  undetermined  constant  C  is  determined  from  empirical  time  domain 
measurements. 


Okal  (1989)  showed  that  the  attenuation  and  dispersion  terms  in  (32)  are  roughly  proportional  as 
follows  when  evaluated  in  the  vicinity  of  20-second  periods'. 


51og(A)  oc  |  log(sin(A))  +  log(e)  ^  log(A) 


(33) 


and 


iog(i/r2)  oc 


(34) 
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where  B  is  a  constant  in  (33).  It  should  be  noted  that  Okal  also  included  a  source  excitation 
correction  in  (34).  However,  (34)  is  a  very  rough  approximation,  which  depends  strongly  on  the 
value  of  dU/dT.  Notice  that  for  Airy  phases,  dU/dT  — *  0,  which  can  cause  large  errors  in  the 
magnitude. 

Substituting  (33)  and  (34)  into  (32)  results  in  the  standard  formulation  for  surface  wave 
magnitudes  (Vanek  et  al,  1962;  von  Seggem,  1977): 

Ms  =  log(a  IT)+B  log(A)  +  C  (35) 

or,  If  only  (34)  is  used  (von  Seggem,  1977;  Rezapour  and  Pearce,  1998): 

Ms  =  log(u  /T)+ -^-log(sin(A))+  BollA + -^log(A)  +  C  (36) 

The  term  Bavin  (36)  is  a  constant  defining  the  attenuation.  It  should  be  noted  that  Rezapour  and 
Pearce  prefer  using  a  coefficient  of  1/3  instead  of  1/2  in  the  dispersion  term  of  (36),  in  order  to 
account  for  distance  effects  of  Airy  phase  propagation. 

5.1  Surface  Wave  Magnitudes  from  Narrow-band  Butterworth  Filters 

The  theoretical  derivation  of  surface  wave  magnitudes  for  narrow-band  filters  is  considerably 
simplified  from  the  above,  since  the  frequency  to  time-domain  transformation  is  non-dispersive, 
for  values  of  coc  satisfying  the  inequality  in  (25),  Recall  the  asymptotic  value  for  an  n*  order 
Butterworth  filter  (equation  23): 

ab=®crnA  (37) 

where  at,  is  the  filtered  time  domain  Butterworth  amplitude,  and  CQt  is  the  Butterworth  comer 
frequency.  Solving  for  A  and  substituting  this  expression  into  (29)  gives  the  corrected 
Butterworth  filtered  time-domain  amplitude  at  period  T: 

rcxA 

abc  =  sin(A)  eUQT  (38) 

tocrn 

Taking  logarithms  of  (38)  results  in  the  time-domain  magnitude  formula: 

Ms(b)  =log(%)+^log(sin(A))+log(e)-^;-log(/c)+Ct  (39) 

where  constants  in  (38)  are  now  combined  in  Cb,  and  coc=  2n  fc.  Notice  that  the  correction  terms 
are  now  geometric  spreading,  attenuation,  and  the  value  of  the  Butterworth  comer  frequency. 
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From  (39)  the  final  formula  for  narrow-band  filtering  with  a  fixed  attenuation  coefficient  can  be 
written  as 


Ms(b)  =  log' K )  +  ^-l°g(sin(A))+  Bau  y-A- log(/c )+  Cb 


(40) 


where  it  is  assumed  that  variations  in  the  attenuation  terms  in  (39)  are  small  around  the  reference 
period  To.  The  constants  Batl  and  Cb  in  (40)  can  be  found  by  empirically  measuring  surface  wave 
amplitudes  a*  across  different  epicentral  distances  A. 


5.2  Normalizing  Butterworth  Magnitudes  to  Standard  Magnitude  Formulae 

An  alternative  to  determining  the  constants  Batt  and  Cb  empirically  is  to  transform  (40)  into  the 
standard  formula  (36)  at  reference  period  To,  and  then  calculate  the  Butterworth  constants  based 
on  standard  formula  constants.  This  has  the  advantage  of  ensuring  that  Butterworth  magnitudes 
are  unbiased  with  respect  to  currently  accepted  magnitudes,  at  least  at  given  reference  periods. 
To  accomplish  this,  first  recall  the  frequency  to  time  transformations  given  in  (22)  and  (23): 


ab  —  6)crnA 


(41) 


The  first  transformation  represents  surface  wave  amplitudes  measured  on  essentially  unfiltered 
and  non-Airy-phase  seismograms,  and  the  second  is  the  amplitude  measured  after  Butterworth 
filtering  with  a  comer  frequency  satisfying  the  inequality  given  in  (25).  Under  these  conditions 
the  unfiltered  time  domain  measurement  can  be  transformed  into  the  filtered  measurement  by 
equating  frequency  domain  amplitudes  in  (41): 


(42) 


Recalling  the  definition  of  a  in  (3),  equation  (42)  can  be  rewritten  in  terms  of  measured 
quantities  as: 


ab  = 


/trV a 


-a 


(43) 


where  G  is  defined  as: 


G  = 


(44) 


k,  A  are  defined  as  in  (29),  and  (oc  =  2nfc.  Also  using  (3),  equation  (25)  can  be  rewritten  as 
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/e* 


G 

t4I 


Notice  that  when  the  equality  holds  in  (45),  ab  =  a  in  (43). 


(45) 


Now,  to  do  the  transformation,  assume  that  a  large  number  of  measurements  have  been  made  at  a 
reference  period  To  on  a  core  seismic  network,  which  is  distributed  across  given  geographic 
areas,  with  source/receiver  propagation  paths  representing  an  average  type  of  geologic  structure. 
Define  the  average  group  velocity  at  To  for  this  structure  as  Uo  and  its  derivative  as  dU/dT\o. 
Define  the  corresponding  value  of  G  in  (44)  as  Go-  Evaluating  (43)  at  Go  and  To  and  substituting 
into  (40)  gives 

Ms  =  log(aro )  +  T-log(sin(A))  +  BattA + ^log(A)  +Cb-  log(G0 )  (46) 

Adding  and  subtracting  log  (To2)  in  (46)  gives 

Ms  =  log(a  /  T0 )  +  ^  log  (sin  ( A)) +  Ban  A  +  ^  log(A)  +  C6  -  log 

\  w  / 

Equating  (36)  and  (47)  shows  that  the  two  equations  are  equal  at  the  reference  period  To  if  C*  is 
defined  as 

C„=C+  log|f-  (48) 

i0 

V. 

and  Batt  is  equivalent  for  both  equations.  Thus,  the  Butterworth  magnitude  (40)  can  be  derived  to 
be  unbiased  with  respect  to  standard  magnitudes  by  defining  Go  and  determining  the  constant  Q, 
from  (48). 


Finally,  it  should  be  noted  that  the  value  of fc  bounded  in  (45)  must  be  calculated  from  actual 
values  of  G,  T,  and  A  for  individual  measurements  in  order  to  minimize  dispersion  error  for  each 
measurement.  However,  since  (45)  is  an  inequality,  if  a  minimum  value  Gmi„  can  be  found  for  all 
expected  propagation  paths,  this  can  be  used  to  define^  as 

f  <  JAhl  (49) 

c  t4 a 

and  this  will  result  in  minimum  dispersion  error  across  the  network  for  all  propagation  paths, 
given  knowledge  of  T  and  A  for  individual  events. 

5.3  Butterworth  Magnitude  Formula  at  Variable  Periods 

It  should  be  emphasized  that  the  above  derivation  is  only  valid  in  the  vicinity  of  20-second 
periods.  The  above  formula  should  not  be  extrapolated  to  short  periods  without  correcting  for 
period-dependent  source  excitation  and  attenuation.  Typically,  explosions  and  shallow 


12 


AFT  AC-TR-04-004 


Theoretical  Analysis  of  Narrow-Band 
Surface  Wave  Magnitudes 


earthquakes  have  source  excitation  functions  which  increase  by  a  factor  of  about  0.2  magnitude 
units  from  20-  to  10-second  periods,  and  again  by  0.2  magnitude  units  between  10-  and  5-second 
periods.  Also,  numerous  studies  (e.g.,  Herrmann  and  Mitchell,  1975)  have  shown  that  the 
attenuation  coefficient  (km'1)  can  increase  by  a  factor  of  three  or  more  from  10-  to  5-second 
periods.  Therefore,  it  is  advisable  to  modify  equation  (40)  with  functions  which  adequately 
approximate  this  behavior,  if  unbiased  short  period  magnitudes  are  required.  However  the 
functions  are  constructed,  they  should  be  normalized  at  T0=  20  seconds  to  reflect  standard 
magnitude  formulae. 

Candidate  functions  can  be  constructed  with  the  form  T0/Tto  normalize  at  T0.  To  modify  the 
attenuation  coefficient,  the  following  form  is  suggested: 


B(T)  =  Bau 


(50) 


When  T  -  To,  the  function  reduces  to  the  constant  coefficient  Batt.  The  value  of  the  coefficient  y 
can  be  chosen  to  reflect  the  increase  the  value  of  the  coefficient  at  shorter  periods. 

For  the  source  excitation,  the  following  function  can  be  constructed  to  correct  the  magnitude 
formula  for  typical  short  period  increases  in  source  amplitude: 


S(T)  =  -S0  log 


(51) 


When  T  -  To,  the  source  function  contributes  no  correction  and,  for  other  periods,  the  value  So 
controls  the  amount  of  correction  at  shorter  periods. 

With  the  above  corrections,  a  final  form  of  the  Butterworth  magnitude  formula  for  variable 
periods  can  be  expressed  as  follows: 


Ms(b)  =  log  (ab  )  + 1  log(sin(A))  +  Ba 


( j1  V 

1  o 

\T  j 


f  T  \ 


A-S0  log 


T 

v 1  J 


-log  (/J+Q 


J'  <  ^min 


tVa 


(52) 

(53) 


To  calculate  MS(b),  the  following  steps  should  be  taken: 

-  Determine  the  epicentral  distance  in  degrees  to  the  event  A  and  the  period  T. 

-  Calculate  the  comer  frequency  fc  of  the  Butterworth  filter  from  (53). 

-  Filter  the  time  series  with  a  zero-phase  Butterworth  band-pass  filter  with  comer  frequencies 

1/T-fc,  1/T+fc 
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ji 


-  Calculate  maximum  amplitude  at  of  filtered  signal.  < 

-  Determine  Ms(b)  from  (52). 

5.4  Specific  Examples 

Three  commonly  used  formulae  for  surface  wave  magnitudes  (amplitudes  measured  zero-to- 
peak,  in  millimicrons)  are  given  by 

Prague  (Vanek  et  al.,  1962): 

Ms  =  log  (a  /  T)  + 1 .66  log(A)  +  0.3  (54) 

Von  Seggem  -  normalized  to  Prague  at  50  degrees  (von  Seggem,  1977): 

Ms  =  log(a  IT) + — log  (sin  (a))  +  0.003 1A  +  |log(A)  +  2.2  (55) 

2  2 

Rezapour  and  Pearce  (Rezapour  and  Pearce,  1998): 

Ms  =  log(a  /  T)  + — log(sin(A))  +  0.0046A  +  ^log(A)  +  2.37  (56) 

2  3 

Plotting  the  magnitude  corrections  as  a  function  of  epicentral  distance  gives: 

Surface  Wave  Magnitude  Corrections 


" "  von  Seggem  (adjusted) 
— '  Prague 


Figure  4. 
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As  can  be  seen  from  the  figure,  the  Rezapour  and  Pearce  correction  is  almost  identical  to  the  von 
Seggem  correction.  Although  the  Prague  formula  biases  magnitude  corrections  to  oceanic 
propagation  paths,  it  is  still  used  to  normalize  newer  formulae  to  maintain  historical  continuity. 
For  purposes  of  this  study,  the  von  Seggem  formula  will  be  used  as  a  baseline,  where  Batt  = 
0.0031,  and  C  =  2.2. 

To  calculate  the  normalized  Butterworth  magnitude  equation,  assume  a  network  defined 
primarily  over  continental  paths,  and  measured  at  the  reference  period  20  seconds.  Let: 

To  =  20  sec 

Uo  =  2.9  km/sec 

(dU/dT)0  =  0.02  km/sec2 

Also  assume  a  3rd  order  Butterworth  filter.  From  equations  (17)  and  (23): 

^ =Xcos(^)=l 

*=1  3 

Using  the  above  values,  and  equations  (44)  and  (48): 


Go  =  0.93,  Cb  =  -0.43 


To  determine  Gmin,  use  equation  (44),  and  look  at  values  of  T,  U,  and  dU/dT  that  minimize  G  for 
various  paths  and  periods  in  the  core  network.  Setting  Gmin=  0.6  should  cover  continental 
signals  between  8  and  40  seconds,  and  oceanic  signals  between  20  and  40  seconds,  including 
mixed  oceanic  and  continental.  This  assumes  a  20-second  oceanic  group  velocity  of  U  = 

3.6  km/sec  and  derivative  dU/dT  =  .08  km/sec2.  The  value  of  Gmi„  should  be  lowered  down  to 
0.2  or  less  for  deep  sediment  structures  at  periods  between  5  and  8  seconds,  and  short  period 
oceanic  paths  between  5  and  20  seconds. 

To  determine  the  period  dependent  attenuation  coefficient,  use  the  form  given  by  (50)  and,  from 
von  Seggem,  use  Batt  =  0.0031  (0.278  x  104  km4).  To  find  y,  find  the  best  fit  of  equation  (50)  to 
empirical  continental  attenuation  coefficient  data.  For  instance,  Herrman  and  Mitchell  (1975) 
published  values  of  anelastic  attenuation  for  the  stable  interior  of  North  America  for  both 
shallow  earthquakes  and  nuclear  explosions.  Plotting  the  values  for  Rayleigh  wave  attenuation 
with  y  =  2.3  gives: 
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ATTENUATION  VALUES  (10**-4  1/km) 


0  10  20  30  40 

Period  (Sec) 

D  D  0  Herrmann  &  Mitchell  (1975) 

Present  Study  0.003 1*(T0/T)**2.3 
1  95%  Confidence  Intervals 


Figure  5. 

Notice  that  the  value  of  B(T)  in  the  present  study  is  lower  than  suggested  by  Herrmann  and 
Mitchell;  this  is  due  to  forcing  the  function  to  fit  Batt  —  0.0031  (0.278  x  10"4  km’1)  at  20  seconds. 
The  point  is  to  ensure  that  B(T)  adequately  reflects  the  significant  increase  in  attenuation  over 
continental  paths  at  short  periods. 

To  determine  the  correction  for  source  excitation,  use  equation  (51)  and  assume  the  source  to  be 
a  shallow  earthquake  or  explosion.  In  this  case,  there  should  be  an  increase  in' amplitude 
between  20  and  10  seconds  by  approximately  0.2  magnitude  units,  with  the  same  increase 
between  10  and  5  seconds.  A  value  of  So  =  0.66  will  correct  for  this  excitation  as  shown  in 
Figure  6: 
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Correction:  -0.66*log(T0/T) 


Figure  6. 

Combining  all  of  the  above  into  (52)  and  (53)  results  in  the  Butterworth  magnitude  formula: 


(20f  \ 

'20^ 

—  A -0.66  log 

T  ) 

T 

\ 1  y 

fc* 


0.6 
rV a 


(57) 

(58) 


Again,  it  should  be  noted  that  a  value  of  Gmi„  =  0.6  should  be  valid  for  most  continental  paths 
with  8  <  T  <  25  seconds,  and  oceanic  paths  with  T  >20  seconds.  For  paths  in  deep  sediments 
with  5  <  T <  8  seconds,  and  oceanic  paths  with  5  <T<20  seconds,  a  value  of  Gmin  =  0.2  should 
be  more  appropriate  to  reflect  the  strong  dispersion  present.  Also  note  that  the  use  of  this 
formula  will  fully  correct  measurement  errors  due  to  Airy  phases. 

6.  Summary 

In  summary,  the  above  theory  gives  a  method  to  calculate  surface  wave  magnitudes  over  broad 
period  and  distance  intervals.  The  fundamental  goals  of  the  study  were  to  develop  a  surface 
wave  methodology  which  can: 


-  Measure  signals  in  the  time  domain  with  minimum  digital  processing  using  Butterworth 
filters. 


-  Effectively  measure  unbiased  surface  wave  magnitudes  at  both  regional  and  teleseismic 
distances,  while  being  applicable  across  a  wide  range  of  periods. 
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-  Ensure  that  the  formula  is  unbiased  with  respect  to  accepted  formulae  at  reference  periods, 
providing  historical  continuity. 

With  rough  estimates  of  group  velocity  and  derivative  bounds,  the  formula  is  unbiased  with 
respect  to  standard  20-second  historical  Ms  measurements,  down  to  5-second  periods  and 
distances  of  less  than  200  km.  The  strength  of  the  method  is  that  it  solves  the  problem  of  trying 
to  simultaneously  measure  in  the  time  domain  well-dispersed  surface  waves  at  teleseismic 
distances  with  non-dispersed  or  Airy  phases  at  regional  distances.  The  weakness  of  the 
methodology  at  this  point  is  the  analysis  of  attenuation.  Although  the  above  formula  (equation 
50)  compensates  for  attenuation  down  to  5  seconds  on  continental  structures,  it  is  clearly 
oversimplified  for  both  regional  and  teleseismic  distances  and  is  the  largest  remaining  source  of 
error  for  surface  wave  magnitudes.  More  study  is  needed  to  define  regionally  varying 
propagation  paths,  including  oceanic,  deep  sediments,  and  basin  and  range  structures,  with 
emphasis  on  analyzing  variability  between  5-25  seconds  across  these  paths.  To  account  for 
regionally  varying  attenuation,  the  theoretical  attenuation  term  in  (39)  can  be  regionalized  as: 


\og(e)XK  ^ 

UQT 


->BV(T)L 


(59) 


where 


B0(T)  = 


\og(e)7TK 

Wwfw 


(60) 


In  (60),  both  U(T)  and  T  can  be  determined  directly  from  individual  narrow-band  filtered 
seismograms;  however,  what  must  be  understood  is  the  variability  of  Q  over  multiple 
source/receiver  paths  ij  at  different  periods  T. 


18 


AFT  AC-TR-04-004 


Theoretical  Analysis  of  Narrow-Band 
Surface  Wave  Magnitudes 


Appendix  A 

Derivation  of  Filtered  Surface  Wave  Integral 


The  purpose  of  this  appendix  is  to  derive  the  solution  to  the  integral 


/„  =-  f  W>  dco 

KJ  0)2n+0)2n 


where 


P  =  t-, 


Start  with  the  1 51  order  Butterworth  filter  integral 


oo 

ri! 


’  .  «c  ei(fico-aco2)  d(Q 

co2  +  co2 


This  integral  can  easily  be  shown  to  be  equivalent  to  the  time  domain  convolution 


2  4i7ta 


which  can  be  expressed  as 


00  XP~r)2 

h=-%=  {e'^e  4a  dT 
24ina  J 


;(P~rf 


iP+tf 


h=  -%=  (e-^e  “a  dr  +  [e-V  J  4a  dr 
2\l  iita  J  J 


The  first  integral  in  (63)  can  be  rewritten  as 


oo  i  2  (iP  V  iP2 
.  — r  -  — +coc  \c+J— 

T  C  4a  l  2a  I  4a  , 
I\  i  =  \  el  v  '  *dz 
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Using  relationships  in  Abromowitz  and  Stegun  (1964),  equation  (64)  has  a  solution 


hi 


uta  eu 


i—+i 
4a 


4amc+—j=: 

24a 


erfc 


4i{4acoc  +  -^L 

2-sfa 


L  V 


(65) 


where  erfc  is  the  complementary  error  function  (1  -  erf). 

Following  the  same  analysis  for  the  second  integral  in  (63)  results  in 


/12  =  -Jin: a  e- 


4£T  [  irfa  J 


4i\  4acor  — i~ 
L  24a 


Adding  (65)  and  (66)  and  substituting  back  into  (63)  gives,  after  rearranging  and  grouping  terms 


/,  =^e  4a 


ef,2  Werf^  (+/?))+ ^erfc^-fi)) 


(67) 


where 


.it 


Vi{±P)  =  e4 


-  :  P  ' 


4act)r±i 


2  4a 


Equation  (67)  is  the  solution  to  the  1st  order  Butterworth  integral  (62).  To  solve  for  the  n*  order 
Butterworth  integral,  expand  out  the  Butterworth  integrand  in  equation  (61)  as  partial  fractions: 


CO, 


2  n 


4 


(0ln  +0)cn  nk=\c°2+4 


(68) 


where 


Mk 


and 


£ k 


7t{2k-n-\) 
2  n 
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Substitute  (68)  into  (61)  for 


f  _£ da) 

ntl’cLe>‘ 


(69) 


The  integral  in  (69)  is  now  equivalent  to  the  1st  order  Butterworth  integral  in  (62)  which  has  the 
solution  (67).  Substituting  (67)  into  equation  (69)  gives  the  final  result: 


4  = 


it 1 

4a  n 


(+^Wc(Tt (+/?)) +eVp*(  P)erfci^fk{-p)) 
znk=\L  ■ 


(70) 


where 


i-t 


'¥k(±/3)  =  e4 


4aek±i-@i= r 

2V« 


—  f°ce 


idk 


a  7t(2k  —  n  —  1) 

“t  = - 

*  In 
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Appendix  B 

Digital  Butterworth  Filters 

Software  for  constructing  digital  Butterworth  filters  is  both  widely  available  and  simple  to 
construct  in  terms  of  recursive  infinite  impulse  response  (HR)  filters.  Kanasewich  (1975)  and 
many  others  discuss  the  construction  of  these  filters  in  detail,  and  it  is  assumed  that  the  reader  is 
familiar  with  these  algorithms.  However,  for  constructing  narrow  band-pass  filters,  it  is 
important  to  point  out  several  issues  in  transforming  continuous  frequency  domain  filters  into 
equivalent  discrete  digital  time  domain  filters  that  can  affect  this  study.  To  accomplish  this,  a 
review  of  the  basic  steps  involved  in  Butterworth  digital  filter  construction  is  appropriate.  * 

B.l  Low-pass  Filter  Design 

a.  Start  with  the  square  of  the  transfer  function  for  a  low-pass  filter  (see  equation  7): 

'  (71) 


It  is  usually  assumed  that  the  square  root  of  (71)  represents  the  amplitude  spectrum  for  a  causal 
one-way  Butterworth  filter.  However,  for  this  study,  only  zero-phase  non-causal  filters  are 
assumed,  realized  by  applying  a  causal  filter  and  then  the  conjugate  reverse  phase  filter,  so 
(71)  represents  the  actual  amplitude  spectrum  for  the  zero-phase  filter. 

b.  Transform  the  filter  into  the  Laplace  domain  with  the  variable  substitution: 


5 


-»  hl  = 


1 


n  In 


l  +  (-l)V 


(72) 


c.  Since  (72)  is  a  real  function,  it  can  be  factored  into  a  polynomial  with  its  conjugate.  The 
polynomial  with  roots  on  the  left  side  of  the  complex  plane  can  be  used  to  construct  a  stable, 
causal  digital  filter: 


hlx  = 


B{s) 


(73) 


d.  To  transform  Hu  into  a  digital  filter,  use  the  bilinear  z  transform  for  a  digital  sampling 
interval  dt : 


dt  1  +  z 

e.  Substitute  (74)  into  (73)  for  an  equivalent  digital  filter 


(74) 
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M 


-  (75) 

k= 1 

Since  z  is  just  the  Laplacian  time  shift  operator,  (75)  can  easily  be  set  up  as  a  recursive  filter  for 
an  input  signal  X(z)  as: 


«  > 

m= 

( »  A 

Xv‘ 

*='  J 

) 

\X(z) 


(76) 


f.  To  run  the  filter,  transform  (76)  into  the  time  domain,  set  initial  conditions,  and  then 
update  y(t)  with  past  values  as  a  recursive  filter.  . 


g.  Finally,  to  run  a  zero-phase  filter,  reverse  the  filtered  output  y(t)  into  x(t),  rerun  (76),  and 
then  reverse  the  filtered  output  y(z)  again.  Again,  notice  that  this  will  have  a  Butterworth  filter 
amplitude  spectrum  given  by  ( 1 ) . 


For  brevity,  the  above  review  is  very  condensed,  and  does  not  include  methods  to  prewarp  the 
Butterworth  cutoff  frequency  coc  to  account  for  non-linearity  in  the  bilinear  z  transform;  and 
methods  to  cascade  the  Laplace  filter  into  smaller  polynomials  to  provide  more  stable  time 
domain  filters.  See  Kanasewieh  (1975)  for  details. 


B.2  Band-pass  Filter  Design 

The  standard  approach  to  constructing  a  band-pass  digital  filter  is  similar  to  the  steps  for  the  low- 
pass  design,  except  for  step  b,  where  a  Laplace  transform  is  substituted  which  maps  a  band-pass 
filter  with  comers  cwi  and  (02  into  a  normalized  low-pass  filter. 

a.  The  transform  that  maps  the  band-pass  into  low-pass  is 

,  00  s  "I-  co^co 2  (77 

c oc  s{cox  -co2) 

b.  Substituting  (77)  into  (71)  gives: 


1 


1  +  (-!)■ 


( s2 

s(^-<o2) 


Y" 

/ 


(78) 


c.  Follow  steps  c  through  g  under  the  low-pass  design  to  realize  the  digital  band-pass  filter. 
Notice  from  (78)  that  when  factoring  the  Laplace  polynomials  into  conjugate  functions  in  (73), 
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there  will  be  a  factor  of  s2n  in  the  numerator  for  the  bandpass  design.  Again,  see  Kanasewich 
(1975)  for  detailed  steps  in  constructing  the  digital  filter. 

B.3  Issues 


One  of  the  problems  in  using  (78)  in  the  construction  of  narrow-band  time  domain  recursive 
filters  is  the  size  of  the  Laplace  polynomial.  From  (78)  it  can  be  seen  that  the  maximum  order  of 
the  polynomial  in  the  denominator  will  be  s4n  instead  of  s2n  as  in  the  low-pass  case.  Although 
some  stability  can  be  gained  by  cascading  filters  down  to  order  2,  numerical  instability  can  be  a 
significant  issue  when  using  single  precision  processing  with  very  narrow  band-pass  filters, 
when  the  filters  are  transformed  into  recursive  form  via  (74)  and  (75).  It  is  essential  that  ’ 
computer  routines  that  design  and  execute  recursive  band-pass  filters  with  a  basic  form  given  by 
(78)  be  double  precision  for  all  floating-point  operators.  In  addition,  for  long  period 
calculations  (>10  sec)  it  may  be  necessary  to  decimate  the  time  series  to  4  samples/sec  or  less  to 
ensure  stability. 

Another  issue  in  using  (78)  is  the  non-linearity  in  transforming  a  symmetric  band-pass  filter  of 
the  form  given  by  (6)  into  the  modified  form  (78).  Recall  that  the  original  band-pass  filter  (6) 
was  transformed  into  a  low-pass  filter  by  a  simple  shift  in  the  frequency  domain  equation  (4). 
Equation  (6)  is 


Hfco)  =  - 


(rn-rn  Y" 


1  + 


CQ-6. ), 

V  ^  , 


and  can  clearly  be  seen  as  symmetric  about  £W0.  Now,  it  can  easily  be  shown  that  (78)  has  the 
following  form  for  its  amplitude  spectrum.  Let  a>\  =  co0  - coc,  and  co2  =  co 0  +  coc.  Then 


H2(co)  = 


1 


1  + 


f  2  \2n 

CO  -CDx(D2 

co(co2-co])^ 


(80) 


It  is  not  obvious,  but  it  can  be  demonstrated  that  H2  is  symmetric  in  log-frequency. 

For  applications  in  this  report,  it  is  instructive  to  plot  (79)  and  (80)  in  a  worst  case  scenario, 
when  the  Butterworth  cutoff  frequency/c  has  the  maximum  width  and  T  is  at  the  shortest  period 
(5  seconds).  Assume  a  3rd  order  Butterworth,  and  let  f0  =  0.2  Hz,/C  =  0.085  Hz.  Plotting  in 
linear  and  log  frequency  gives: 
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Amplitude  Spectrum  -  Linear  Frequency 


e  0.2  0.4  0.6 

Frequency  (Hz) 

—  HI 
H2 


Amplitude  Spectrum  -  Log  Frequency 


Frequency  (Hz) 

—  HI 
"  H2 


Figure  B-l. 


From  Figure  B-l ,  it  would  appear  that  there  is  a  significant  effect  due  to  how  the  band-pass  filter 
is  constructed;  however,  this  is  only  pronounced  due  to  expanding  small  changes  with  the 
logarithmic  scale  on  the  amplitude  axis.  Fast  Fourier  transforming  the  above  frequency  filters 
into  the  time  domain  to  give  the  impulse  response  show  that  for  purposes  of  measuring 
magnitudes,  the  above  effects  cause  only  2n“  order  changes  in  the  time  domain: 


Impulse  Response 


Time  (Sec) 

—  HI 
--  H2 


Figure  B-2. 

Again,  the  above  is  the  worst-case  scenario.  As  the  period  increases  and  the  filter  width 
decreases,  it  can  be  shown  that  the  difference  between  H\  and  Hi  decreases.  In  addition,  running 
synthetic  seismograms  with  both  source  excitation  and  attenuation  built  in  show  less  than  .01 
magnitude  difference  when  filtered  with  H\  and  Hi. 
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B.4  Alternative  Time  Domain  Recursive  Filter 


As  an  alternative  to  the  standard  HR  design  given  above,  a  different  transformation  can  be  used 
to  design  a  recursive  narrow-band  filter  which  has  the  following  advantages  for  operational 
processing: 

-  The  filter  response  is  equivalent  to  (79),  following  the  exact  theory  in  this  report. 

-  The  filter  is  built  as  cascaded  first  order  filters  instead  of  quadratic  filters,  thus  improving 
the  precision  for  very  narrow-band  applications. 

-  The  filter  returns  a  complex  time  domain  response,  which  can  be  used  to  extract  the  real 
filtered  response  and  its  envelope  function  without  resorting  to  Fast  Fourier  transforms. 

-  The  algorithm  has  been  successfully  tested  for  narrow-band  applications  in  this  report,  for 
periods  between  5  and  40  seconds,  and  sampling  rates  up  to  40  samples/second,  without 
requiring  decimation. 

The  basic  steps  to  construct  the  filter  are: 

a.  Transform  the  low-pass  into  band-pass  as: 


.  a  s-io0 

p  =  i—  = - - 

0)c 


(81) 


or  equivalently 


s  =  pa)c+ico0  (82) 

This  formula  represents  a  simple  translation  of  a  low-pass  filter  from  the  origin  to  a>0,  resulting  in 
a  complex  band-pass  filter.  Note  that  this  is  similar  to  using  the  positive  half  of  the  Fourier 
spectrum  to  form  a  complex  time  series,  in  order  to  calculate  the  Hilbert  transform  for  the  time 
series  envelope  function  (Papoulis,  page  124,  1962). 

b.  Substitute  (81)  into  (7 1 )  for 


Hb  -> 


i+HrO’f 


/ 


1 +  H)" 


S  —  iCOn 


\2" 


\  / 


(83) 


c.  Equation  (83)  can  be  factored  into  poles  by  first  finding  the  poles  pj  for  the  prototype 
low-pass  filter  in  (83)  (Kanasewich,  page  181,  1975)  and  then  substituting  these  poles  into 
equivalent  band-pass  poles  in  (82),  resulting  in  the  form 
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q)2/s 

Hb=- — s. - ,  Sj  =PjO)c+iC0o  (84) 

n^-5y) 

d.  Equation  (84)  can  then  be  cascaded  into  simple  first  order  filters 

H„  =-^~  (85) 

S-Sj 

which  can  be  transformed  into  time  domain  recursive  filters  of  the  form  (76)  using  the  bilinear 
z-transform  (74). 

e.  Finally,  after  successively  running  the  cascade  recursive  filters  and  then  reversing  for 
zero  phase,  the  real  part  of  the  complex  time  series  can  be  extracted  for  the  narrow-band  filtered 
signal,  and  the  modulus  of  the  series  can  be  calculated  for  the  envelope  function. 

B.5  Alternative  Recursive  Filter  Coding  Algorithm 

To  realize  the  above  filter,  a  simple  algorithm  can  be  constructed  which  can  be  easily  coded  into 
Cor  FORTRAN  code: 

-  Define  the  following  variables  and  vectors  •  _  _ 

INTEGER: 

m  =  Butterworth  order 
n  =  number  of  time  series  points 
j,k  =  counters 

REAL 

dt  =  sampling  interval 

wo  =  Butterworth  band-pass  center  frequency  {2itfo ) 

coc  -  Butterworth  band-pass  comer  frequency  {2rfc) 

xr(n)  —  input  unfiltered  time  series 

yr(n)  =  output  filtered  time  series 

er(n)  =  output  envelope  time  series 

ermx  =  maximum  value  of  envelope  series 

COMPLEX  DOUBLE  PRECISION 

p(m)  -  prototype  low-pass  Butterworth  poles 
s(m)  =  equivalent  band-pass  Butterworth  poles 
a\(m)  =  Z-transform  recursive  coefficients 
al(m)  =  Z-transform  recursive  coefficients 
z\(n)  =  complex  time  series 

z2(n)  =  complex  time  series 
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Calculate  complex  poles  of  Butterworth  polynomial 


Pj  =  exp  —~(2j  —  \  +  m) 
_2m 


j  =  1,2...  m 


sj  =Pj(»c+icoQ 

-  Calculate  complex  bilinear  z-transform  recursive  coefficients 


ah=-^L- 

J  2 -sdt 


j  =  \,2...m 


2  +  s  ,dt 

a2.  = - 2— 

J  2-s.dt 


-  Place  input  time  series  xr  into  real  part  ofzl 


zlk  =  xrk 


k  =  1,2... n 


-  Calculate  m  cascaded  first  order  filters 


j  =  1,2 ...m  ; 

z2k  =  zh 
zl,  =alyz2, 


k  =  l,2...n 


z^k  alj(z2k  +z2k-i)+a2jzh-i  k  —  2,3...n 
-  Reverse  complex  time  series 


Z2  k  —  Z^n-k+ 1 


k  =  1,2  ...n 


Calculate  m  reversed  first  order  filters  (al*,a2*  complex  conjugates  of  al,a2) 


j  =  l,2...m  ; 
zh  =  z2k 
z 2,  =  ul*  zl, 


k  =  1,2  ...n 


z2  k  (zl*  +  zl*_i )+ <z2^.z2i_,  k  —  2,3  ...n 
-  Reverse  complex  time  series 


Z^k  ~  z2„-k+\ 


k  =  l,2...n 
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-  Calculate  output  real  time  series,  envelope,  and  envelope  maximum 

yrk  =  2  REAL(zlk)  k  =  \,2..M 

erk=2\zlk\ 

ermx  =  MAX(erk  ,ermx ) 

NOTE:  For  applications  given  in  this  report,  bilinear  z-transform  pre-warping  (Kanasewieh, 
page  192,  1979)  is  not  required for  the  above  recursive  algorithm,  due  to  the  low-pass  filter 
being  designed  well  below  the  Nyquist  frequency. 
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